ChIP AML PiPeline v2¶
In [1]:
import os
import pandas as pd
import sys
sys.path.insert(0, '../..')
import itertools
from scipy import stats
import numpy as np
from JKBio.epigenetics import chipseq as chip
from JKBio.utils import helper
from JKBio.utils import plot
import igv
import SimpSOM as sps
from scipy import stats
import seaborn as sns
from matplotlib import cm
from matplotlib import pyplot as plt
from IPython.display import IFrame
import seaborn as sns
from bokeh.plotting import *
import igv
import numba
from numba import jit
from scipy.cluster.hierarchy import linkage, leaves_list
from sklearn.cluster import AgglomerativeClustering, DBSCAN, KMeans, OPTICS
from sklearn.mixture import GaussianMixture
from sklearn.manifold import MDS, TSNE
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from IPython.display import IFrame
from statsmodels.stats.multitest import multipletests
from pybedtools import BedTool
import pyBigWig
output_notebook()
%load_ext autoreload
%autoreload 2
In [21]:
#setting basic parameters
project="Cobinding_ChIP"
version="v3"
merging_version="motifs"
window="150"
In [32]:
#retrieving the data from the previous notebook
%store -r mergedmot
%store -r cols
%store -r annot
%store -r version
%store -r window
%store -r crc
merging_version
Out[32]:
Correlation with Annotations¶
non binarized¶
In [35]:
for val in mergedmot.columns[72:]:
mergedmot.loc[mergedmot[mergedmot[val]==0].index,val.split('_')[0]]=0
merged = mergedmot[mergedmot.columns[:72]].rename(columns={val.split('_')[0]: val.split('_')[0]+"_filtered" for val in mergedmot.columns[72:]})
In [41]:
#data = pd.DataFrame(np.corrcoef(stats.zscore(0.01+merged[merged.columns[cols:]].T)), columns=merged.columns[cols:], index = merged.columns[cols:])
data = pd.DataFrame(np.corrcoef(merged[merged.columns[cols:]].T.astype(bool)), columns=merged.columns[cols:], index = merged.columns[cols:])
link = linkage(data[:annot-cols])
col = data.iloc[annot-cols:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
rdb = cm.get_cmap('RdBu_r', 256)
for val in col.columns:
a = [rdb(128+int(v*128)) for v in col[val]]
col[val] = a
In [42]:
#correlation_withannotation
fig = sns.clustermap(data.iloc[:annot][data.columns[leaves_list(link)]], vmin=-1,vmax=1, col_colors=col.T, figsize=(10,20), cmap='RdBu_r', row_linkage=link, col_cluster=False, colors_ratio=0.01)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_withannotation")
#fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_correlation_with_annotation.pdf')
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_boolcorrelation_with_annotation.pdf')
plt.show()
binarized¶
In [43]:
#data = pd.DataFrame(np.corrcoef(stats.zscore(0.01+merged[merged.columns[cols:]].T)), columns=merged.columns[cols:], index = merged.columns[cols:])
data = pd.DataFrame(np.corrcoef(merged[merged.columns[cols:]].T.astype(bool)), columns=merged.columns[cols:], index = merged.columns[cols:])
link = linkage(data[:annot-cols])
col = data.iloc[annot-cols:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
rdb = cm.get_cmap('RdBu_r', 256)
for val in col.columns:
a = [rdb(128+int(v*128)) for v in col[val]]
col[val] = a
In [44]:
#correlation_withannotation
fig = sns.clustermap(data.iloc[:annot][data.columns[leaves_list(link)]], vmin=-1,vmax=1, col_colors=col.T, figsize=(10,15), cmap='RdBu_r', row_linkage=link, col_cluster=False, colors_ratio=0.008)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_onbinarized_signal_withannotation")
#fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_correlation_with_annotation.pdf')
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_boolcorrelation_onbinarized_signal_with_annotation.pdf')
plt.show()
Looking at peak overlap¶
How many of peak in A (column) overlaps with peak in B (rows)
in other words:
what is the percentage of B's peaks that are overlaped by A's peaks
In [45]:
overlap, correlation, pvals, enrichment = chip.pairwiseOverlap(merged, norm=True)
In [46]:
#saving the enrichments
overlap.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'pairwise_overlap.csv')
correlation.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'pairwise_correlation_on_overlap.csv')
enrichment.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'pairwise_enrichment_on_overlap.csv')
pvals.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'pairwise_enrichment_on_overlap_pvals.csv')
pairwise overlap¶
In [47]:
data = pd.DataFrame(data=overlap,index=merged.columns[cols:], columns=merged.columns[cols:])
link = linkage(data.iloc[:annot-cols]) # D being the measurement
c = data[[co for co in data.columns if co not in data.iloc[:annot-cols].index.tolist()]]
viridis = cm.get_cmap('viridis', 256)
for val in c.columns:
a = [viridis(int(v*255)) for v in c[val]]
c[val] = a
In [48]:
fig = sns.clustermap(data.iloc[:annot-cols][data.columns[np.concatenate((leaves_list(link), range(annot,len(data.columns))))]], row_linkage=link, col_colors=c,figsize=(20,18), colors_ratio=0.01, col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("pairwise_overlap_clustermap")
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'pairwise_overlap_clustermap.pdf')
plt.show()
Correlation on overlaps¶
on the overlaps given above
In [49]:
data = pd.DataFrame(data = correlation, index = merged.columns[cols:], columns = merged.columns[cols:])
link = linkage(data.iloc[:annot-cols], optimal_ordering=True) # D being the measurement
c = data[[co for co in data.columns if co not in data.iloc[:annot-cols].index.tolist()]]
viridis = cm.get_cmap('viridis', 256)
for val in c.columns:
a = [rdb(128+int(v*128)) for v in c[val]]
c[val] = a
In [50]:
fig = sns.clustermap(data.iloc[:annot-cols][data.columns[np.concatenate((leaves_list(link), range(annot,len(data.columns))))]], row_linkage=link, col_colors=c, colors_ratio=0.01, figsize=(20,18), vmin=-1, vmax=1, cmap='RdBu_r', col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_on_pairwise_overlap_clustermap")
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_correlation_onoverlap.pdf')
plt.show()
Enrichment¶
In [51]:
link = linkage(enrichment, optimal_ordering=True)
fig = sns.clustermap(enrichment,figsize=(20,20), vmin=-5,vmax=5, cmap='RdBu_r', col_linkage=link, row_linkage=link)
fig.ax_col_dendrogram.set_visible(False)
fig.ax_row_dendrogram.set_visible(False)
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_enrichment_clustermap_all_to_all.pdf')
plt.show()
In [52]:
enrichment.loc[enrichment.columns[leaves_list(link)]][enrichment.columns[leaves_list(link)]].values.max()
Out[52]:
In [53]:
link = linkage(enrichment, optimal_ordering=True)
plot.correlationMatrix(enrichment.loc[enrichment.columns[leaves_list(link)]][enrichment.columns[leaves_list(link)]].values, dataIsCorr=True, names=enrichment.columns[leaves_list(link)].tolist(), pvals=pvals[enrichment.columns[leaves_list(link)]].loc[enrichment.columns[leaves_list(link)]].values, size=12, title= version + '_' + merging_version + '_' + window + 'enrichment', folder='../results/'+project+'/plots/',interactive=True,minval=-4,maxval=4)
Out[53]:
looking only at crcs¶
In [54]:
enrichment = enrichment[crc].loc[crc]
In [ ]:
link = linkage(enrichment, optimal_ordering=True)
plot.correlationMatrix(enrichment.loc[enrichment.columns[leaves_list(link)]][enrichment.columns[leaves_list(link)]].values, dataIsCorr=True, names=enrichment.columns[leaves_list(link)].tolist(), pvals=pvals[enrichment.columns[leaves_list(link)]].loc[enrichment.columns[leaves_list(link)]].values, size=12, title= version + '_' + merging_version + '_' + window + 'enrichment', folder='../results/'+project+'/plots/',interactive=True,minval=-4,maxval=4)
clustering¶
I have tried gaussian mixtures and Agglomerative clustering algorithm. Only the second can create a hierarchical clustering.
It seems that gaussian mixture makes more sense given the data we have, for now, is more "homogeneous".
I am still not so happy with the clustering. It can be because of the how much importance, outlier values and the high number of noisy values from locations with no peaks.
We can use similar methods to RNAseq to improve this (clamping values, log transform, first round of PCA..)
using KMeans¶
In [96]:
data= merged[merged.columns[cols:annot]].values
# using score
scaled_data = stats.zscore(data)
# using regular scaling
#scaled_data = (data-data.min(0))/(data.max(0)-data.min(0))
In [56]:
data.max(0)
Out[56]:
In [97]:
#https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html#sklearn.cluster.KMeans
n = 1
n_clust=[10,20,50,100,300]
kmean = KMeans(n_clusters=n_clust[n],n_jobs=8)
groups = kmean.fit_predict(scaled_data)
centroid = kmean.cluster_centers_
plotting center location of clusters (kind of enrichment)¶
In [116]:
Out[116]:
In [111]:
df/
Out[111]:
In [109]:
df = pd.DataFrame(centroid,columns=merged.columns[cols:annot]).T
In [118]:
df = pd.DataFrame(centroid,columns=merged.columns[cols:annot]).T
#df = (df.T/df.max(1)).T
fig = sns.clustermap(df, vmax=1)
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_kmeans_'+str(n_clust[n])+'_centroids.pdf')
showing the cobinding matrix with the clusters¶
In [61]:
merged = merged.reset_index(drop=True)
In [119]:
rand = np.random.choice(merged.index,5000)
subgroups = groups[rand]
sorting = np.argsort(subgroups)
rainb = cm.get_cmap('gist_rainbow', len(set(groups)))
colors = [rainb(i) for i in subgroups[sorting]]
viridis = cm.get_cmap('viridis', 256)
data = merged[merged.columns[annot:]]
for val in data.columns:
data[val]=np.log2(1+data[val])
data = (data-data.min())/(data.max()-data.min())
data = data.iloc[rand].iloc[sorting]
for val in data.columns:
a = [viridis(int(v*256)) for v in data[val]]
data[val] = a
data["clusters"] = colors
In [120]:
fig = sns.clustermap(np.log2(1.01+merged[merged.columns[cols:annot]].iloc[rand].iloc[sorting].T), vmin=0, vmax=3, figsize=(20,20), z_score=0, colors_ratio=0.01, col_cluster=False,col_colors=data, xticklabels=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("sorted clustermap of cobindings clustered with Kmeans")
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_kmeans_'+str(n_clust[n])+'_clustermap_cobinding.pdf')
plt.show()
Computing enrichment on the clusters¶
In [121]:
enr, pvals, _ = chip.enrichment(merged, groups=groups)
In [122]:
# enrichment for each cluster
fig = sns.clustermap(enr.T,figsize=(14,20), col_cluster=0, vmax=4, vmin=-4, cmap='RdBu_r')
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_kmeans_'+str(n_clust[n])+'enrichment_on_cluster_metasubset_zcorescale.pdf')
plt.show()
using bokeh¶
pvalue will be informed by opacity
In [67]:
fig = plot.correlationMatrix(enr.loc[enr.columns[leaves_list(link)]][enr.columns[leaves_list(link)]].values, maxokpval=10**-9, pvals= pvals[enr.columns[leaves_list(link)]].loc[enr.columns[leaves_list(link)]], names=enr.columns[leaves_list(link)].tolist(), size=12, title= version + '_' + merging_version + '_' + window + 'enrichment', folder='../results/'+project+'/plots/', interactive=True, minval=-4, maxval=4, dataIsCorr=True)
showing enrichment over cobinding¶
In [68]:
plot.andrew(groups, merged, annot, cols=8, precise=False, folder = '../results/' + project + '/plots/' + version + '_' + merging_version + '_' + window + '_kmeans_', title = "sorted clustermap of cobindings clustered with Kmeans")
Using SOMs¶
In [69]:
#Import the library
size = 20
#Build a network 20x20 with a weights format taken from the raw_data and activate Periodic Boundary Conditions.
net = sps.somNet(size,size, scaled_data, PBC=True)
#Train the network for 10000 epochs and with initial learning rate of 0.01.
net.train(0.01, 10000)
#Save the weights to file
net.save('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'_cobinding_SOMweights_'+str(size))
In [712]:
net = sps.somNet(0,0, scaled_data, loadFile='../results/'+project+'/'+version+'_'+merging_version+'_'+window+'_cobinding_SOMweights_'+str(size), PBC=True)
In [70]:
pathplot = "../results/"+project+'/plots/'+version+'_'+merging_version+'_'+window+"_somNets/"
! mkdir $pathplot
In [86]:
#Print a map of the network nodes and colour them according to the first feature (column number 0) of the dataset
#and then according to the distance between each node and its neighbours.
for col in range(scaled_data.shape[1]):
net.nodes_graph(colnum=col, printout=True, show=True, path="../results/"+project+'/plots/'+version+'_'+merging_version+'_'+window+"_somNets/", colname=data.columns[col])
In [87]:
diffs = net.diff_graph(show=False, returns=True)
plt.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_cobinding_SOM_'+str(size)+'.pdf')
In [89]:
plot.SOMPlot(net, size=size, colnames=merged.columns[cols:annot], folder='../results/' + project + '/plots/' + version + '_' + merging_version + '_'+window, minweight=0.05)
In [90]:
#Cluster the datapoints according to the Quality Threshold algorithm.
clusts = net.cluster(scaled_data, type='KMeans',numcl=n_clust[n])
sorting = []
for clust in clusts:
sorting.extend(clust)
sorting = np.array(sorting)
In [91]:
viridis = cm.get_cmap('gist_rainbow', len(clusts))
colors=[]
for i, clust in enumerate(clusts):
colors.extend([viridis(i)]*len(clust))
colors = np.array(colors)
viridis = cm.get_cmap('viridis', 256)
data = merged[merged.columns[annot:]]
for val in data.columns:
data[val] = np.log2(1+data[val])
data = (data-data.min())/(data.max()-data.min())
rand = np.random.choice(range(len(sorting)),5000)
rand.sort()
data = data.iloc[sorting[rand]]
for val in data.columns:
a = [viridis(int(v*256)) for v in data[val]]
data[val] = a
data["clusters"] = [tuple(i) for i in colors[rand]]
In [92]:
#sorted clustermap of cobindings clustered with SOM
fig = sns.clustermap(np.log2(1.01+merged[merged.columns[cols:annot]].iloc[sorting[rand]].T), vmin=0,vmax=3,figsize=(20,15),z_score=0, colors_ratio=0.01, col_cluster=False,col_colors=data, xticklabels=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("sorted clustermap of cobindings clustered with SOM")
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_SOM_'+str(n_clust)+'_clustermap_cobinding.pdf')
plt.show()
In [93]:
#let's look at the enrichment over SOM clusters
groups = []
for i, val in enumerate(clusts):
groups.extend([i]*len(val))
enr, pvals, _ = chip.enrichment(merged, groups = np.array(groups))
In [95]:
fig = sns.clustermap(enr.T,figsize=(12,12), vmax=4, vmin=-4, cmap='RdBu_r')
fig.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_kmeans_'+str(n_clust[n])+'enrichment_on_SOMcluster_metasubset.pdf')
plt.show()
Plot TSNE density map¶
In [123]:
rand = np.random.choice(merged.index,30000)
In [124]:
red_data = TSNE(2,600,verbose=10,n_iter=1500).fit_transform(scaled_data[rand])
np.save('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'red_data.npy',red_data)
In [125]:
# Let's look at the data. the density map of the distribution of all pseudo enhancer over their tf binding profile
_, ax = plt.subplots(figsize=(10,10))
fig = sns.kdeplot(red_data[:,0], red_data[:,1], shade=True, ax=ax)
plt.savefig('../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+'_density_TSNE.pdf')
In [126]:
# now let's look at it from anothe rploting method. we can see many many islands
plot.bigScatter(red_data,binsize=0.4,showpoint=False,precomputed=False, logscale=True, title='density plot of enhancers in TF cobinding space with TSNE', folder='../results/'+project+'/plots/'+version+'_'+merging_version+'_'+window+"_")
In [ ]:
## TODO: color by dependencies/pathways/expression.. after ABC model,